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ABSTRACT 


The  computer  program  for  the  solution  of  unsteady  heat  conduction 
problems  Involving  plane  or  axisymmetric  geometry,  devised  by  Robert  E. 
Nickell,  has  been  modified  to  include  temperature  dependent  properties 
and  time  dependent  boundary  conditions.  The  original  IBM  7094  computer 
dependent  program  has  been  converted  for  use  on  the  IBM/OS  360  Model  67 
computer.  In  both  programs  FORTRAN  IV  language  vas  used.  A  "User's 
Manual"  has  been  constructed  for  the  modified  program. 
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I.  INTRODUCTION 


l 

|  „ 


*» 


In  this  chapter  the  reasons  and  advantages  of  the  finite  element 
analysis  method  are  given.  The  purpose  of  this  presentation  and  the.  ob¬ 
jectives  of  the  study  are  discussed. 

A.  REASONS  FOR  THE  FINITE  ELEMENT  ANALYSIS 

With  the  ever  increasing  complexity  of  problems  occuring  from  engi¬ 
neering  situations  in  design  and  analysis,  there  is  a  rising  need  for  an 
approximate  method  of  solution  for  these  problems.  This  technique!  must 
be  flexible  enough  to  encompass  a  large  variety  of  problems  and  st..ll  give 
accurate  answers  to  the  simplest  of  situations.  Sines  most  approximate 
methods  of  analysis  of  complex  engineering  problems  Involve  many  tedious 
calculations,  it  would  enhance  the  usefulness  of  the  technique  if  It  were 
capable  of  being  programed  for  computer  solution.  This  would  relieve 
the  engineer  of  the  tedious  burden  of  hand  calculated  solutions  and  leave 
him  more  time  to  devote  to  the  job  of  engineering. 

The  finite  element  method  of  analysis  is  an  approximate  method  of 
analysing  engineering  problems  which  meet  these  specifications.  It  has 
a  high  degree  of  flexibility  and  can  be  readily  applied  to  computer  sol¬ 
utions  to  obtain  rapid  solutions  to  complex  problems. 

B.  OBJECTIVES  OF  THIS  STUDY 

The  major  objective  of  this  study  is  a  computer  program  using  the 
finite  element  method  of  analysis,  for  the  solution  of  two-dimensional 
unsteady  heat  conduction  problems.  In  fulfilling  this  goal  there  were 
several  minor  objectives  accomplished  along  the  way.  The  original  pro¬ 
gram  was  written  for  the  IBM  7094  system.  It  was  converted  for  operation 
on  the  IBM/OS  360  computer  at  the  Naval  Postgraduate  School  Computer  Center. 
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Involved  In  this  conversion  was  the  changing  of  control  cards  and  a  few 
FORTRAN  IV  language  differences.  Also  faster  random  access  disks  were 
substituted  for  the  magnetic  tapes  originally  used  for  temporary  storage. 
The  input  and  output  formats  were  overhauled  to  make  better  use  of  space 
on  the  printed  output  page.  More  extensive  labels  were  added  to  the 
printed  output  in  order  to  ease  the  use  of  the  program  by  persons  not 
familiar  with  the  finite  element  technique. 

C,  PURPOSE  OF  THE  PRESENTATION 

The  purpose  of  this  presentation  is  to  make  available  a  computer 
program  using  the  finite  element  method  of  analysis  for  the  solution  of 
unsteady  heat  conduction  problems.  In  doing  so  an  introduction  into  the 
formulation  of  the  finite  element  method  is  given  for  those  not  familiar 
with  the  technique,  to  give  an  indication  of  the  power  of  this  type  of 
analysis.  A  description  of  the  components  of  the  program  is  included  to 
enable  the  user  to  better  understand  the  logic  oc  the  computer  solution 
and  increase  the  user's  ability  to  formulate  problems  to  which  the  pro* 
gram  can  be  applied.  Specific  instructions  on  formulation  of  problems 
and  loading  of  input  data  ia  made  available  in  a  chapter  on  user  infor¬ 
mation.  Finally,  as  la  customary  in  works  of  this  type,  recommendations 
for  future  modifications  of  the  program  and  possible  applications  are 
Included. 
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II.  THE  FINITE  EUSMENT  FORMULATION 

In  tnis  chapter  the  finite  element  method  is  formulated.  The  idea 
of  a  functional  of  a  function  and  the  variation  of  a  functional  ore  Intern* 
duced.  Theae  principle*  are  developed  into  an  approximate  eolation  of 
the  two-dimensional  unsteady  heat  conduction  equation  using  e  matrix 
iterative  solution. 

A.  FORMULATION  OF  THE  VARIATIONAL  METHOD 

If  a  functional  of  a  temperature  field  can  he  found  such  that  Its 
extremum  yields  an  equation  describing  unsteady  heat 


fc'jIjt=ecT-i 


<1-0 


and  an  equation  describing  the  heat  flux  across  the  surface  of  the  * 


n<  fcoTij  =  H* 


<a-D 


then  the  variation  of  the  functional  will  yield  a  set  of  first  order, 
ordinary  differential  equations  in  terms  of  the  nodal  temperatures.  DsUg 


the  method  of  tensor  notation,  a  comma  denotes  differentiation  with  res¬ 
pect  to  the  following  subscript  and  repeated  subscripts  imply  aisesstliiii. 

Let  IT  be  a  functional  of  the  temperature  field  T(x,y,*,t)  and  the 
first  time  derivative  of  the  temperature  field  T(x,y,s,t),  he  defined  hy 


fct,  fcg  n*cTt-  iT)d/  ♦  UiVT4s 


<2-J) 


The  *  variation  of  7T(T,T)  with  respect  to  T  {with  T  held 


la  glvan  by 


<Stt^  <Stt(t  W,T) 
<$€ 


(1-4) 


13 


where  £  is  a  small  parameter  and  A  is  one  of  a  family  of  functions. 

A  is  zero  on  the  boundary  of  the  region  and  arbitrary  elsewhere.  An 
extremum  of  the  functional  Tr(T,T)  implies  that  £tt(t  ,  T)  equals  zero,  i 


<ftrtr+frA,T) 

it 


6  -O 


O 


(2-5) 


Fir^t 


TKT+ (.-!,+)  = 


(iCr+tA^fujCIW)^ 


0 

fc.(J n^T+t^dS 

\)o 


(2-6) 


Using  the  identity 


i  |(T*^),tK?j(r+«),JJ=  Jfr*  a),u  -fttj(T-*fcA),jL 


A 


if  the  conductivity  tensor  is  symmetric  i.e.  kij  «*  kj£.  Then 

Cr a,t)  =  x 

KtAT'^A)dv  - 


(2-7) 


The  volume  Integral 


Integral. 


j(T,;  ktjA^v  can  be  transformed  into  a  surface 


(T,;K;jA),j  4v  »  ^  nitty 


(2-8) 


.e. 
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When  (2-7)  is  evaluated  at  £  *»  0 


)  -  J (cru. A)rj  -  Kij  i>jc^ 

+ec4t-i/C)<dv  AA3 

<5tt(t^-6A,t)  =  j  (-k^  +ect~i)AJv 

f  *£^1  * > ;  Xds  -  \  n£aA<te 
d$  J  J3 

The  extremum  of  TT  (^reduce*  to 

^Cj%jC  = 


in  the  region  V.  Alao  on  the  surface  S 


(2-*> 


<2-10) 


«^yT-.j  -  nin 


a-u) 


These  two  equations  (2-10,  2-11)  are  the  unsteady  heat  conduction  oqnatia 
and  the  boundary  ;lux  equation.  Thia  verifies  that  a  function  T  which 
yields  an  extremum  of  the  functional  defined  by  (2-3)  satisfies  both  the 
unsteady  heat  conduction  equation  In  the  region  V  and  the  boundary  flux 
equation  on  the  surface  S  QQ*. 

B.  FORMULATION  OF  THE  FINITE  ELEMENT  PROBLEM 

» 

If  the  choice  of  T  and  T  is  restricted  so  that  they  are  represented 
by  certain  constants  which  are  obtained  in  their  formulation,  the  func¬ 
tional  TT  becomes  a  real-valued  function.  In  this  caa*"TT  shell  be 

where  and  ft  are  vectors  of  nodal  peint  value* 

of  The  extremum  of  must  no*  be  found  which 


numbers  in  brackets  refer  to  similarly  numbered  items  in  Bibli¬ 


ography. 
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implies 


=■  O  (2-12) 

cSTf 

From  this  point  on  the  finite  element  formulation  will  be  considered 
in  two  dimensions  only.  The  region  V  will  be  divided  into  a  number  plane 

m 

triangular  elements.  The  functions  T  and  T  can  be  uniquely  described 
within  each  of  these  elements  by  the  values  of  the  function  at  each  of 
the  nodal  points  of  the  element  T^,  Tj  »  Tm»  and  linear  function  of  the 
co-ordinates  of  the  element. 

Let  the  temperature  field  in  an  element  be  given  by 


TU,v,t)  = 


(2-13) 


and  the  time  rate  of  temperature  change  given  by 

t^.V,t)  =  <«S(x,v»[A][t?(t)]  (2.14) 

is  a  vector  which  specifies  the  special  approximation  of  T  and  T. 
The  matrix  [V]  i*  *  constant  and  is  defined  by  the  above  relationships. 
Its  exact  value  and  derivation  will  be  discussed  in  detail  in  section  II, 
E. 

The  temperature  gradient  T,^  ( fj  ■  x,y)  can  be  expressed  in  terms  of 
the  nodal  temperatures  by 


(2-15) 


In  the  following  development  the  subscript  f  ,  denoting  nodal  point 
values,  will  be  drepped  and  will  b 2  assumed. 
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Writing  the  functional  IT  in  terms  of  ncdal  point  values. 


-  WWVi)«»v  ~  1 

0 


(2-16) 


Taking  the  variation  of  TT  with  respect  to  ^T^  and  setting  It  equal  to 
sero 


^tt(  w,  m) = [*>}  -(f} +[a©-  &} =o  (2-17) 


where 


w*  l  bWwmmjv 

ov 

[c]  -  ^  pcMV<d>M«Jv 

[i*]A  iWVdv 

Jv 

{%*}  =[  WV)%|JS 


(2-18) 


(2 -m 


(2-20) 


(2-21) 


C.  BOUNDARY  CONDITIONS 

There  are  four  types  of  boundary  conditions  to  be  considered.  (1) 
specified  temperature,  (2)  specified  flux,  (3)  convection  boundary  layer, 
and  (4)  adiabatic.  For  a  specified  temperature  boundary  condition  is 
constant  over  boundary  segment  Sj_.  For  a  specified  flux  q  »  TJ  over 
boundary  segment  S2.  For  a  convection  boundary  layer  q  -  h(Ti  -  T^)  over 
boundary  segment  S3.  Finally  for  an  adiabatic  boundary  condition  q  ••  0 
over  boundary  segment  S4, 
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The  boundary  integral  (2-21)  becomes 


where 


{%*}  -in  *  Mm  -  ri 


{ %'}  - 


-X  /  ,^T 


(2-22) 

(2-23) 


[H]  -U  [A]V<<2S>MdS 


(2-24) 


1  KXl  [A]W<)S  (2-25) 

The  Integral  on  segment  is  sero  because  the  variation  of  the  functional 
was  specified  as  sero  on  the  boundary.  The  integral  on  segment  S4  is 
sero  since  there  is  no  heat  flow  across  the  boundary. 

To  obtain  the  extremum  of  the  functional  IT  the  nodal  point  tempera¬ 
tures  ^T]  which  satisfy  the  following  matrix  equation  must  be  found  [d}. 

(M-M)St]  +  [c]{t]  -  [%*}  °  (2-20 


D,  THE  MATRIX  EQUATION 

Let  the  time  variable  be  ti  such  that  t^  »  i&k  i 
For  this  development  i  will  refer  to  cime  increments. 


0,1,2,... 


Let 


M  •  KHh) 


(2-27) 


Equation  (2-17)  can  be  written  as 

W[t}-  -{fV, 


(2-28) 
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vhere 


(fk-rMri-M 


(2-29) 


Assuming  that  is  constant  over  the  Interval  t^  £  t  t^  +  j,,  t.a,  r 

IT  h  -  (ft*  ~lrk)M 

lT};. ,  can  be  obtained  by  a  Taylor's  expansion  around  t  »  t^. 


tT}i+l  =iTh  f^k 


(2-30) 


Hit.,  >iTU  mifh  +  ^([rii+, -Itk)  (2-3i) 


Mm  ={T}t  *^(lrL,  ♦  {T}.J 


(2-32) 


Using  equations  (2-29)  and  (2-32)  ««*  can  be  detsreiaad 

In  terms  of  and  {l}.  .  Solving  (2-32)  for  gives 


Vrh,  -  aCFlu  -ltk)  -  {% 


(2-33) 


Substituting  this  Into  (2-28)  gives 


H  +  E^Airii  *(%)  <«* 


Also  from  equation  (2-28) 


mi  =  celt!  - 


(2-35) 
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Substituting  (2-34)  into  (2-35)  gives 


-if)  vl4., 


(2-36) 


or 


(2-37) 


In  a  multi-element  body  the  element  matrix  equations  must  be  assembled 
into  a  single  matrix  equation.  The  assembly  and  solution  of  the  equation 
is  a  complex  task  but  it  can  be  greatly  simplified  by  the  use  of  a  com¬ 
puter  [jSH . 


E.  FORMULATION  OF  THE  ELEMENT  MATRIX 

If  the  region  to  be  analysed  is  divided  into  plane  triangular  elements 
then  a  typical  element  vould  appear  like  the  following  figure. 


The  temperature  in  an  element  can  be  described  in  terms  of  a  linear  com¬ 
bination  of  the  coordinates  of  the  elements  nodal  points  and  the  valves 
of  the  temperature  at  each  of  the  nodal  points. 

Assuming  that  the  temperature  distribution  is  of  the  form 

T(X,VL)  =  +0(i(±.)X  +<=<3(.-Oy  (2-38) 
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The  temperature  at  each  of  the  nodal  points  can  be  written  as 


Ti.  -  °T  + 

Tj  ~  °'1  +  xj°<2-  *•  Vj°<3  (2-39) 

Xrv  =  V 

Solving  for  the  coefficients  2’  ant*  °^3  using  Cramer**  rule  the 

following  results  are  obtained. 


I  Xi  Yt 

L&t-  ,  Yi  - 

i  x*  V 

Note:  is  the  area  of  the  triangular  element. 


(2-40) 


<*,  = 

^5  “ 


Ti 

yj 

y± 

T» 

/ 

Tt 

1 

Ti 

Xi 

i 

T*. 

Y* 

1 

Ki 

J 

Ai 

1 

•m 

IA 

IA 

ZA 
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°i, = (Tteb  t  t  v?  i + 1^ 

^41;  ih:  iM.'  i)/* 

V(h|,  j+  h\\  4  TJl  Jjl)/ZA 


Recalling  that  the  temperature  distribution  was  previously  described  by 


the  expression 


= 


(2-42) 


is  the  vector  of  the  nodal  temperatures  and  [a]  is  the  spacial 


dependence  matrix  which  is  a  constant.  Let 


<£>  =  <i  k  v> 


(2-43) 


It}- 


4,2.  4,3 
4z3 

°lj<  °( J2-  °<33 


=  <0>[a]It} 
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80 


c\ 


(%  U  +  9fc)j  +  VW 

T  -  <(/  X  7)  J(%Ti  +  -v 

« 

(*>»,  1c  +  Vli+  *33*1 


J 


T-  fat,  l^-f-  %,T./j  T*)  +  •f°feTDx  +t%3trJj  4  pfjjTi^'y 

Therefore,  comparing  (2-43a)  with  (2-41) 
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The  spacia?  dependence  matrix  is  then  defined  by 
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F.  THE  METHOD  OF  LIMPING  ELEMENT  PROPERTIES 

When  working  with  a  computer  it  is  usually  more  convenient  to  work 
with  a  quadrilateral  element.  The  quadrilateral  is  divided  into  four 
triangles  with  the  center  nodal  point  common  to  each  of  the  triangles. 
(See  figure  2-2). 


Figure  2-2  Quadrilateral  Element 
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The  matrix  equation  (2-28)  can  be  assembled  for  this  element  in  the 
following  form 
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Partitioning  this  equation  results  in 
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This  equation  can  be  multiplied  out  to  yield  two  equations 


(2-45) 


+  T5  +  \\\  +  ?5  •  (£4  <*-“> 

<KSJ>  iT4  +  K55T5  +  +  °55?s  *  f5  «’47> 

The  specific  heat  matrix  JjC  J  is  approximated  by  lumping  the  heat  ca¬ 
pacities  at  the  four  external  nodes.  Thus  |j[Q  becomes  a  diagonal  matrix 
(4x4)  and  C55  a  0.  Therefore  (2-46)  can  be  written  as 


M  h]  +  t5  +  Mi  it}  ■  [£l]  «-M> 

and  (2-47)  as 

<>'5J>  {Ti}  +  %5t5  -  h  «-*’) 

Solving  (2-49)  for  T^  and  substituting  into  (2-48)  results  in 

(fill  W  +  (Kis)  K55  (£5  -  <K5j>[T?j)  +  [Cij](*l}  -  f£i]  (2*50) 
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or 


(  IX?  -  i\il  K55  <hj)  >  {Ti}  +  tiJXi  ■  ffi]  -  Xis}  %5f: 


Equation  (2-51)  is  analogous  to  equation  (2-28)  except  that  [V]  and 
are  now  (4x4)  matrices  and  is  a  (4x4)  vector  L0- 
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III.  PROGRAM  STRUCTURE 


This  chapter  describes  the  general  composition  of  the  program  and' 
gives  a  description  of  the  function  of  each  section  of  the  program. 

A.  PROGRAM  IDENTIFICATION 

The  program  was  written  by  Robert  E.  Nickell  while  in  the  Department,  of 
Civil  Engineering,  University  of  California  at  Berkeley.  It  wae  formerly 
identified  as  W2498,W035~2— 57343,  Nickell.  It  was  written  in  July  1968 
and  was  revised  by  Allan  B.  Chaloupka  in  May  1969.  Formulation  was 
based  upon  a  variational  principle  derived  by  M.  E.  Gurtin  of  Brown 
University  QQ. 

The  program  has  been  renamed  DFETHC-Finite  Element  Transient  Heat 
Conduction. 

B.  PURPOSE 

The  purpose  of  the  program  is  to  produce  high  speed  solutions  to 
transient  heat  conduction  problems.  These  problems  may  involve  plane-  or- 
axi symmetric  geometry,  non-homogeneous  anisotropic  material  configura¬ 
tions,  temperature  dependent  material  properties,  and  time  dependent 
boundary  conditions. 

C.  PROGRAMING  INFORMATION 

The  program  was  originally  written  in  FORTRAN  IV  for  the  IBM  7094 
computer.  It  was  converted  for  use  on  the  IBM  OS/360  Model  67  computer. 

The  program  was  also  rewritten  in  double  precision  (REAL  *8). 

D.  PROGRAM  CAPACITY 

The  sise  problem  that  can  be  analysed  by  this  program  is  subject  to 
the  following  limitationa" 

Maximum  number  of  nodal  points  . . ,  „ .  500 
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Maximum  number  of  quadrilateral  elements  .  490 

Maximum  number  of  materials  .  25 

Maximum  difference  of  nodal  point  numbers 

for  the  same  e  lement  . . .  31 


These  are  not  stringent  limitations  and  they  may  be  adjusted  some¬ 
what  by  the  user  to  fit  individual  problems.  The  limits  are  changed  by 
changing  the  appropriate  DIMENSION  and  statements  and  one  card  in 

the  segment  that  calculates  the  half -band  width.  The  only  real  limita¬ 
tion  placed  upon  the  sice  of  the  problem  that  can  be  analysed  is  the 
amount  of  core  storage  available  in  the  computer  used. 

E.  GENERAL  STRUCTURE 

The  program  consists  of  a  main  program  and  several  subprograms  used 
for  repetitive  operations.  The  main  program  can  be  divided  into  three 
parts.  The  first  part  is  concerned  with  reading  into  the  program  and 
printing  out  the  data  describing  each  problem.  It  is  portrayed  by  the 
flow  diagram  in  figure  3-1. 

The  second  portion  of  the  main  program  is  concerned  with  the  routing 
of  operations  according  to  certain  logical  variables.  These  variables 
describe  four  possible  combinations.  They  are  constant  material  proper¬ 
ties  and  constant  boundary  conditions,  temperature  dependent  material 
properties  and  constant  boundary  conditions,  constant  material  proper¬ 
ties  and  time  dependent  boundary  conditions,  and  temperature  dependent 
material  properties  and  time  dependent  boundary  conditions.  Each  of  these 
possibilities  is  portrayed  by  a  flow  diagram  in  figures  3-2  through  3-5. 

The  third  portion  of  the  main  program  is  concerned  with  calculating 
the  effective  forcing  vector,  solving  the  matrix  equation  and  printing 
the  temperatures  for  each  time  increment.  The  operation  is  then  directed 
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Figure  3-1  Main  Program  Flow  Diagram,  Part  I 
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Figure  3-2  Main  Program  Flow  Diagram,  Part  II 
Constant  Properties  and  Boundary  Conditions 
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Figure  3-3  Main  Program  Flow  Diagram,  Fart  II 
Temperature  Dependant  Propertiaa  and 
Constant  Boundary  Conditions 


31 


Figure  3-4  Main  Program  Flow  Diagram,  Part  II 
Constant  Properties  and  Time  Dependent  Boundary  Conditions 


32 


Figure  3-5  Main  Program  Flow  Diagram,  Part  II 
Temperature  Dependent  Properties  and 
Time  Dependent  Boundary  Conditions 
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back  to  the  beginning  of  the  second  section  of  the  main  program  to  start 
calculations  for  the  next  time  increment.  The  third  portion  of  the  main 
program  is  portrayed  by  a  flow  diagram  tn  figure  3-6. 

F.  MAIN  PROGRAM 

The  main  program  is  composed  ot  a  number  of  individual  steps  which 
comprise  the  finite  element  solution  technique,  and  a  system  for  the  in¬ 
put  and  output  of  problem  information.  Not  included  in  the  main  program 
are  the  matrix  assembly  processes,  the  matrix  equation  solving  processes, 
and  the  boundary  condition  input  routines,  Theie  are  included  in  sub¬ 
programs  FORM, SYMSOL, TFUNC, CFUNC, FCBC, FTBC,  and  FFBC. 

1.  Control  Information 

The  first  information  to  be  read  into  the  program  is  the  control 
information.  This  information  sets  the  basic  variables  for  each  of  the 
steps  of  the  program  and  also  routes  the  problem  to  various  sections  of 
the  program.  This  information  includes  the  number  of  nodal  points,  the 
number  of  elements,  the  number  of  convection  boundary  conditions,  the  num¬ 
ber  of  materials,  the  number  of  time  sequences,  whether  the  problem  is  in 
plane  or  axisymmetric  geometry,  the  output  print  interval,  the  system  of 
units,  whether  a  final  spacial  temperature  distribution  is  desired,  and 
the  reference  temperature  of  the  body. 

Also  read  in  with  the  control  information  are  two  titles.  The 
first  is  a  72  space  title  that  can  be  used  by  the  user  to  label  each 
problem.  The  second  is  a  48  space  title  that  is  used  to  input  the  labels 
used  with  the  system  of  units  chosen. 

The  appropriate  titles  and  control  information  are  printed  on  the 
first  page  of  the  printed  output,  for  each  problem. 
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2.  Material  Properties 

The  thermal  conductivity,  density,  specific  heat,  and  heat  gener¬ 
ation  properties  for  each  material  are  read  into  the  program.  The  materi¬ 
al  properties  are  then  printed  in  order  on  the  printed  output  page  and 
labeled  with  the  appropriate  units.  Since  the  program  is  for  anisotropic 
mater}  jls,  kxx»kyy>  an<*  ^xy  are  rea<^  in  order  to  form  the  thermal  con¬ 
ductivity  tensor.  The  specific  heat  and  density  are  read  in  as  a  product 
in  order  to  reduce  the  number  of  variables  in  the  program. 

3.  Time  Sequencing  Information 

In  this  section  the  time  sequencing  information  is  read  in  for 
each  sequencing  scheme.  This  information  consists  of  two  logical  varia¬ 
bles  and  two  numerical  variables.  The  logical  variables  indicate  whether 
there  are  temperature  dependent  properties  and  if  there  are  time  dependent 
boundary  conditions.  The  two  remaining  variables  describe  the  number 
of  time  steps  and  the  sice  of  each  time  increment.  The  time  sequenci  g 
information  is  printed  on  the  printed  output  page  with  appropriate  labels. 

4.  Nodal  Point  Information 

Ordinarily  this  segment  of  the  program  will  read  the  x,y  coordinates 
and  initial  temperature  of  each  nodal  point  and  print  the  information 
with  the  appropriate  dimensions.  Because  of  the  possibility  of  a  large 
number  of  nodal  points  (500)  and  the  equally  large  possibility  of  human 
error  in  punching  data  cards  there  is  a  segment  of  the  program  that  will 
generate  nodal  point  information.  This  segment  has  a  few  limitations, 
which  will  be  discussed  in  chapter  IV,  but  it  is  usually  able  to  generate 
almost  all  of  the  interior  nodal  points.  As  an  example  in  the  configura¬ 
tion  shown  in  figure  3-7  the  nodal  points  indicated  by  squares  were  read 
into  the  program  and  the  remaining  nodal  points  were  generated  by  the 
program. 
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nodal  Point  ar.d  Element  Generation 


5.  Element  Information 


The  purpose  of  this  segment  is  about  the  same  as  the  previous 
one.  The  number  of  the  element,  the  number  of  each  of  its  nodes,  the 
type  of  material,  and  the  amount  of  heat  generation  within  the  element 
are  read  into  the  program  and  printed  out  with  appropriate  labels.  Again 
this  segment  has  the  capability  to  facilitate  the  loading  of  data  by 
generating  element  information.  This  segment  also  has  some  limitations 
which  will  be  discussed  in  chapter  IV,  but  it  is  usually  able  to  generate 
most  of  the  elements  in  a  given  configuration.  As  an  example,  in  figure 
3-7  all  the  elements  were  generated  except  those  that  are  starred. 

Also  included  in  this  segment  is  a  routine  that  calculates  the 
half -band  width  of  the  problem.  This  is  an  important  control  parameter 
that  is  used  later  in  the  program  to  assemble  the  matrix  equation. 

6.  Time  Steps 

The  time  step  portion  of  the  program  includes  two  large  loops. 

The  first  loop  executes  the  program  through  each  time  sequence  of  the 
problem.  The  second  loop  executes  the  program  through  the  required  num¬ 
ber  of  time  Increments  for  each  time  sequence.  Inside  these  two  loops 
the  two  logical  time  sequencing  variables  and  the  loop  counters  control 
the  routing  of  the  problem. 

The  two  logical  variables  control  the  choice  of  variable  material 
properties  and  boundary  conditions.  If  the  material  properties  and  bound¬ 
ary  conditions  are  constant  then  the  coefficient  matrix  and  forcing  vector 
are  calculated  only  once.  The  matrix  equation  is  then  solved  for  suc¬ 
cessive  time  increments  using  only  the  original  matrices.  If  variable 
material  properties  are  involved  the  effective  coefficient  matrix  must 
be  reapplied  since  the  coefficient  matrix  is  reformed  from  zero  values. 

If  time  dependent  boundary  conditions  are  involved  then  the  contribution 
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of  the  boundary  conditions  to  the  coefficient  matrix  must  be  reapplied. 
The  coefficient  matrix  from  the  original  material  properties  is  .kept  in 
temporary  storage.  When  both  time  dependent  boundary  conditions  and 
temperature  dependent  properties  are  used  then  the  coefficient  matrix 
and  forcing  vectors  must  be  reformed  for  each  time  incraattnt  based  upon 
the  ianediate  time  and  temperature  values  and  Initial  properties.  The 
data  for  temperature  dependent  properties  is  calculated  from  the  initial 
conditions  and  the  present  nodal  temperatures.  The  data  for  time  depend - 
ent  boundary  conditions  is  either  read  into  the  program  with  -each  time 
increment  or  generated  within  the  program. 

Different  time  sequences  can  be  run  for  the  same  problem.  This 
enables  the  user  to  speed  or  slow  the  time  change  during  a  particular 
aegamnt  of  the  problem. 

7.  Cyprectigq  .Bgmyla^  Cgnditijpa 

A  convection  boundary  condition  occurs  when  the  -problem  includes 
a  convection  boundary  layer  as  part  of  its  boundary  condition*.  One  con¬ 
vection  boundary  condition  consists  of  two  nodal  point  numbers,  the  ‘Con¬ 
vective  heat  transfer  coefficient  and  the  aablent  fluid  temperature.  The 
two  nodal  points  define  the  convection  boundary  for  each  condition. 

Time  varying  convection  boundary  conditions  can  be  obtained  by 
varying  the  adbient  fluid  temperature.  The  new  fluid  temperature  may  be 
reed  with  each  new  time  step  or  it  may  be  generated  by  wing  the  function 
subprogram  FCBC.  FCBC  is  used  in  the  same  manner  ss  the  function  sub¬ 
programs  for  temperature  dependent  material  properties.  The  argument  of 
FCBC  ie  TIME  end  the  value  of  FCBC  is  currently  set  at  one.  If  tempera 
ture  dependent  material  properties  are  used  then  the  convection  boundary 
condition  section  stores  and  resppliss  the  convection  boundary  condition 
contributions  to  the  coefficient  matrix  and  the  forcing  vector. 
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8.  Temperature  and  Flux  Boundary  Conditions 

A  temperature  or  flux  boundary  condition  occurs  when  the  problem 
specifies  a  temperature  or  flux  at  a  nodal  point.  One  temperature  or  flux 
boundary  conditions  consists  of  the  nodal  point  number  and  the  value  of 
the  temperature  or  flux  at  that  nodal  point. 

Time  varying  temperature  and  flux  boundary  conditions  can  be  ob¬ 
tained  by  varying  the  temperature  or  flux  at  each  nodal  point.  The  new 
temperature  or  flux  may  be  read  with  each  new  time  step  or  it  may  be 
generated  by  using  the  function  subprograms  FTBC  or  FFBC  for  time  depend¬ 
ent  temperature  or  flux  respectively.  FTBC  and  FFBC  are  used  in  the 
same  manner  as  the  function  subprograms  for  temperature  dependent  materi¬ 
al  properties.  The  argument  of  each  is  TIME.  Their  current  values  are 
both  one.  If  temperature  dependent  properties  are  used  the  temperature 
and  flux  boundary  condition  section  stores  and  reapplies  the  boundary 
conditions  to  the  coefficient  matrix  and  the  forcing  vector. 

9.  Temperature  Output 

This  segment  of  the  program  is  concerned  with  the  print  out  of 
the  temperature  solutions  for  each  nodal  point  at  each  time  step.  The 
value  of  the  time  is  printed  and  the  nodal  point  numbers  and  correspond¬ 
ing  temperatures  for  that  time  are  printed  in  rows  below.  If  it  is  de¬ 
sired,  the  temperature  can  be  printed  with  special  data  at  the  end  of 
each  time  sequence. 

G.  SUBROUTINE  FORM 

The  subroutine  FORM  is  used  to  calculate  the  effective  element  ther¬ 
mal  conductivity,  heat  capacity,  and  heat  generation.  The  subroutine 
performs  these  operations  for  each  of  the  elements  of  the  configuration 
and,  if  temperature  dependent  properties  are  us  sd  j  for  cod'  t iiSv  incr Smcut  • 
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Figure  3-8  The  Matrix  Equation 


Figure  3-9  Computer  Storage  for  the 
Matrix  of  Coefficients 
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First  the  mean  x,y  coordinates  and  the  mean  temperature  of  the  eleamnt 
are  calculated  by  averaging  the  values  of  each  of  the  nodes.  Using  the 
Man  x,y  coordinates  the  quadrilateral  element  is  divided  into  four  tri- 
angler.  The  conductivity  tensor,  heat  capacity,  and  heat  generation 
vectors  are  formed  for  each  of  the  four  triangular  sections.  As  these 
matrices  are  formed  they  are  assembled  into  the  larger  (5)  nodal  point 
equation  for  the  quadrilateral  element.  The  heat  capacity  and  heat  gen¬ 
eration  vectors  are  assembled  so  that  the  fifth  term  of  each  vector  is 
dropped  when  forming  the  lumped  properties  matrix  equation.  The  center 
nodal  point  is  eliminated  from  the  thermal  conductivity  tensor  (5x5), 
according  to  equation  (2-52),  to  give  the  thermal  conductivity  tensor 
^4x4)  for  the  lumped  properties  matrix  equation.  The  thermal  conductivity 
tensor  and  heat  capacity  and  heat  generation  vectors  for  each  elaamnt 
are  returned  to  the  main  program  where  they  are  loaded  into  the  coeffi¬ 
cient  matrix  and  vectors  used  in  the  solution  of  the  total  matrix  equa¬ 
tion  for  each  time  step.  For  problems  where  the  properties  of  the 
materials  are  not  temperature  dependent  the  thermal  conductivity  tensor 
and  the  heat  capacity  and  heat  generation  vectors  are  calculated  only  at 
the  beginning  of  the  problem.  If  temperature  dependent  material  proper¬ 
ties  are  included  these  matrices  are  calculated  at  the  beginning  of  each 
time  step. 

H.  SUBROUTINE  SYMSOL 

The  subroutine  SYMSOL  in  conjunction  with  the  loading  amthods  of  the 
main  program  is  used  to  solve  the  finite  eleamnt  matrix  equation.  It 
consists  of  two  parts.  The  first  part  reduces  the  coefficient  matrix 
to  its  final  form.  The  second  part  forms  the  effective  forcing  vector 
and  solves  for  the  temperatures  at  each  of  the  nodal  points. 
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In  an  ordinary  formulation  the  final  matrix  equation  for  the  finite 
element  method  would  look  like  figure  3-8.  The  matrix  of  coefficients 
g|  is  a  banded  symmetric  matrix.  In  the  computer,  in  order  to  conserve 
core  storage,  the  matrix  of  coefficients  is  stored  in  columns.  The  first 
column  is  the  diagonal  of  the  matrix  of  coefficients.  (See  figure  3-9). 
All  the  operations  performed  on  the  matrix  of  coefficients  are  done  on 
the  upper  half  of  the  symmetric  banded  matrix  in  its  stored  form. 

After  the  coefficient  matrix  is  formed  it  is  put  into  upper  triangular 
form.  This  yields  a  matrix  equation  of  the  form  shown  below. 


Figure  3-10  Matrix  Equation  in  Upper  Triangular  Form 

This  is  accomplished  by  the  first  segment  of  SYMSOL,  The  second  portion 
of  SYMSOL,  using  the  effective  forcing  vector,  starts  with  the  temperature 
at  the  bottom  of  the  nodal  temperature  vector  and  back  substitutes  to 
solve  for  the  nodal  temperatures.  These  operations  are  performed  for  each 
time  Increment. 
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I.  TEMPERATURE  DEPENDENT  PROPERTIES 

There  are  two  function  subprograms  that  are  used  in  conjunction  with 
temperature  dependent  properties.  They  are  TFUNC  and  CFUNC.  TFUNC  is 
used  for  the  temperature  dependent  thermal  conductivity  and  CFUNC  is  used 
for  the  temperature  dependent  product  of  the  density  and  specific  heat. 

In  the  programs  present  configuration  each  of  these  function  subprograms 
is  prepared  for  temperature  independent  properties  and  t'n*.ir  values  are 
therefore  equal  to  one.  If  the  thermal  conductivity  were  temperature 
dependent  in  the  form  oft  k  *  kQ(l  +  bT),  where  b  is  a  constant  coef f i- 
'ient  and  k0  is  the  value  of  the  thermal  conductivity  at  T  .  0,  then 
TFUNC  would  have  to  be  modified  by  the  user  to  include  the  factor  1  +  bT, 
The  material  properties  read  into  the  program  would  then  be  aero  tempera¬ 
ture  properties. 
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IV.  USER  INFORMATION 


This  chapter  describes  the  manner  in  which  problems  should  be  formu¬ 
lated.  Also  the  procedures  used  in  inputting  data  into  the  program  are 
inc luded . 

A.  INPUT  DATA 

Data  for  each  problem  processed  by  this  program  is  assembled  in  groups. 
The  first  group  includes  the  titles  and  all  the  control  parameters.  The 
second  group  is  the  material  properties.  The  third  group  is  the  time  se¬ 
quencing  information.  The  fourth  group  is  nodal  point  information.  The 
fifth  group  is  element  information.  The  sixth  group  is  boundary  condi¬ 
tions.  After  each  data  deck  there  must  be  three  blank  cards. 

B.  CONTROL  PARAMETERS 

The  following  control  parameters  are  read  into  the  program. 

NUMNP  number  of  nodal  point 

NUMEL  number  of  elements 

NUMCBC  number  of  convection  boundary  conditions 

NUMMAT  number  of  materials 

NSEQ  number  of  time  sequences 

KAT  type  of  geometry,  0  implies  plane  geometry 

INTER  output  print  interval 

UNIT  units  system 

KPUNCH  special  temperature  output 

JUMP  generated  temperature  of  flux  boundary  conditions 

KEY  generated  convection  boundary  conditions 

TO  initial  temperature 
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Also  included  with  the  control  data  are  two  labels.  TITLE  (18)  is  an 
arbitrary  label  that  may  be  used  to  identify  individual  problems.  It  is 
72  spaces  long.  The  second  label  is  UN(12).  UN(12)  is  an  array  of  four 

unit  labels.  Each  has  12  spaces.  The  units  must  be  printed  in  the  fol¬ 
lowing  order:  length,  mass,  time,  and  temperature.  These  labels  must  be 
left  justified. 

C.  MATERIAL  PROPERTIES 

The  following  material  properties  are  read  into  the  program. 

XCOND  k 

xx 

YCOND  kyy 

XYCOND  k 

xy 

SPHT  £>c 

HX  heat  generation 

Material  properties  are  read  in  with  one  material  to  a  card.  If  the 
material  is  isotropic  k^  -  kyy  -  0.  If  the  material  is  anisotropic  the 
three  values  will  rot  be  the  same.  Each  material  is  given  a  number  ac¬ 
cording  to  the  order  it  is  read  into  the  program. 

If  temperature  dependent  material  properties  are  used  both  logical 
variables  must  be  true.  Also  both  JUMP  and  KEY  must  be  some  number  other 
than  zero  and  FCBC,FTBC,  and  FFBC  must  be  one. 

D.  TIME  SEQUENCING  INFORMATION 

The  following  time  sequencing  information  is  read  into  the  program. 

LTAG  temperature  dependent  properties  (true  or  false) 

Ml AG  time  dependent  boundary  conditions  (true  or  false) 

ITAG  number  of  time  steps 

TAG  time  increment 

The  time  sequencing  information  controls  the  operation. of  the  program 
through  the  use  of  the  logical  variables  LTAG  and  MTAG  but  it  can  also 
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be  used  to  slow  or  speed  time  during  the  problem.  An  example  would  be 
two  sequences  with  the  first  solving  for  the  temperature  60  times  during 
the  first  hour  of  the  problem  (i.e.,  every  minute).  The  second  time  se¬ 
quence  would  solve  for  the  temperature  6  times  in  the  second  hour  of  prob¬ 
lem  »;itue  (i-  e. ,  every  10  minutes) . 

E.  NODAL  POINT  CONSTRUCTION 

The  following  nodal  point  Information  is  read  into  the  program  for 
nodal  point  data. 

N  number  of  the  nodal  point 

KODE  1  indicates  a  temperature  boundary  condition 

X  x  coordinate  of  the  nodal  point 

Y  y  coordinate  of  the  nodal  point 

T  initial  temperature  of  the  nodal  point 

The  data  for  each  nodel  point  is  put  on  one  card. 

Earlier  it  was  mentioned  that  the  program  could  generate  nodal  points. 
This  feature  has  a  few  limitations  that  must  be  observed  or  incorrect 
information  will  result.  A  nodal  point  whose  KODE  is  1  cannot  be  gener¬ 
ated.  Nodal  points  must  be  generated  a  Ling  s',  reight  lines.  Tie  initial 
temperature  of  the  generated  nodal  points  must  be  ine  same  and  equal  to  TO, 
the  initial  temperature.  To  generate  nodal  points  simpi>  . sve  r  the 
data  cards  between  two  nodal  points  on  a  straight  line  and  the  program 
will  generate  the  information  for  the  nodal  points  in  between. 

In  constructing  the  mesh  of  nodal  points  for  a  problem  it  is  sometimes 
impossible  to  divide  the  area  to  be  analysed  into  uniform  rectangular 
quadrilaterals.  At  times  it  will  not  be  desirable  to  do  so.  Consider 
figure  4-1.  In  preparing  a  nodal  mesh  for  this  configuration  the  first 
step  would  be  to  assume  that  by  symmetry  there  is  no  heat  flux  across  the 
center  lines.  This  reduces  the  figure  to  one  of  its  quadrants  as  !:• 
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Figure  4-2  Reduced  Irregularly  Bounded  Configuration 


figure  4-2„  The  else  of  the  quadrilateral  elements  would  vary  depending 
upon  hear  interested  the  user  was  in  a  certain  tret  of  the  figure.  In 
figure  4-2  the  area  around  the  curve  la  being  aore  closely  examined  but 
the  aesh  else  aust  also  be  snail  to  enable  the  user  to  approacinate  any 
curved  surface  using  a  series  of  straight  line  segaants.  It  is  almost 
inpossible  to  divide  an  area  with  a  curved  boundary  without  encountering 
a  right  triangular  element.  To  convert  this  to  a  quadrilateral  elenent 
place  a  nodal  point  in  the  sriddle  of  the  hypotenuse. 


Figure  4-3  Conversion  of  a  Right -Triable 
to  a  Quadrilateral 

In  mistering  the  nodal  points  it  aust  be  kept  in  mind  that  the  aaai- 
atuzc  difference  in  nodal  point  numbers  in  one  elenent  is  31.  This  Is  a 
limitation  set  by  the  half-band  width  of  the  computer  program.  It  Is 
convenient  to  Increase  the  numbering  of  the  nodal  points  and  the  eloasnt* 
in  the  same  pattern,  although  it  is  not  a  requiresient. 

F.  ELEMENT  CONSTRUCTION 

The  following  elenent  information  la  read  into  the  program  for  ele¬ 
ment  data. 

N  number  of  the  element 

1X(N,K)  K  -  1  -  4  numbers  of  the  element**  nodal  points 
EC(N,5)  material  type 

HTGEN  element  heat  generation 

The  date  for  each  element  must  be  placed  on  one  card.  The  ntdal  point 
numbers  are  placed  in  clockwise  order  starting  with  the  nuafeer  in  the 
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upper  left  hand  corner.  Element*  may  be  generated  by  the  program  but 
there  are  some  limitation*.  Only  element*  of  the  same  material  type  can 
be  generated  at  one  time.  The  heat  generation  must  he  the  same  for  the 
elements  generated  together.  Elements  must  be  generated  in  rows  or  col¬ 
umns.  To  generate  elements  simply  leave  out  those  element  cards  that 
you  wish  generated. 

G.  BOUNDARY  CONDITIONS 

There  are  two  types  of  boundary  conditions  that  may  be  read  into  the 
program,  convection  and  temperature  or  flux.  (An  adiabatic  boundary  con¬ 
dition  is  a  flux  boundary  condition  with  the  flux  equal  to  aero.)  The 
following  information  is  read  into  the  program  for  convection  boundary 
conditions. 

IB  nodal  point  number 

JB  nodal  point  number 

HB  convective  heat  transfer  coefficient 

TEMPB  ambient  fluid  temperature 

The  data  for  each  portion  of  the  convection  boundary  is  put  one  one  card. 
If  time  varying  ambient  temperature  is  to  be  generated  by  the  program  the 
control  parameter  KEY  must  be  some  number  other  than  zero,  say  one.  If 
the  time  varying  ambient  temperature  is  not  generated  by  the  program  new 
data  for  each  segment  of  the  convective  boundary  must  be  read  into  the 
program  with  each  new  time  step.  In  this  case  KEY  is  zero. 

The  following  information  is  read  into  the  program  for  temperature 
and  flux  boundary  conditions. 

NN  nodal  point  number 

TQ  temperature  or  flux 

The  data  for  the  temperature  or  flux  boundary  condition  at  each  nodal 
point  is  put  on  one  card.  If  tine  varying  temperature  or  flux  is  to  be 
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generated  by  the  program  the  control  parameter 'JUMP  must  be  some  number 
other  than  zero.  If  time  varying  temperature  or  flux  is  not  generated 
by  the  program  new  data  for  each  boundary  nodal  point  muat  be  read  Into 
the  program  with  each  new  time  step.  In  this  case  JIMP  is  sera.  Both 
temperature  and  flux  boundary  conditions  must  be  introduced  into  the  pro¬ 
gram  in  the  same  manner.  For  both  convection  and  temperature  or  flat 
boundary  conditions,  if  time  variations  are  used  the  logical  variable 
MIAS  must  be  true.  Also,  regardless  of  the  position  of  the  last  nodal 
point,  a  boundary  condition  card  must  be  read  into  the  program  for  it. 

If  the  last  nodal  point  is  not  a  boundary  point  then  the  data  card  con¬ 
sists  of  the  nodal  point  number  and  a  value  of  zero  in  the  temperature 
or  flux  field. 

H.  UNITS 

The  units  used  for  each  problem  must  be  consistent  with  the  input 
data  or  the  labels  on  the  output  will  be  incorrect.  Presently  there  are 
six  different  units  systems  available.  Each  system  consists  of  a  length, 
mass,  time,  and  temperature.  Each  system  has  a  number  which  must  be  read 
in  as  a  control  parameter  (UNIT).  The  systems  are: 

1.  feet,  pound  mass,  hours,  °F 

2.  inches,  pound  mass,  seconds,  °F 

3.  feet,  pound  mass,  seconds,  °F 

4.  meters,  kilograms,  seconds,  °C 

5.  centimeters,  grams,  seconds,  °C 

6.  feet,  pound  mass,  minutes,  °F 

If  the  user  does  not  prefer  the  units  available  he  can  use  another 
system  in  his  input  and  change  the  labels  or  just  ignore  the  labels  al¬ 
together.  The  system  of  units  will  depend  greatly  upon  the  configuration 
of  the  problem  and  the  temperature  differences  involved. 
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I.  ERROR  EXITS 


There  are  three  error  exits  in  the  program.  The  first  occurs  if  a 
non-positive  number  of  nodal  points  is  read  into  the  program.  This  exit 
is  used  to  stop  the  program.  The  other  two  occur  in  the  nodal  point  and 
element  information  reading  sections.  If  these  statements  are  executed 
then  the  execution  of  the  program  is  terminated.  If  either  of  the  last 
two  exits  are  executed  a  statement  will  be  printed  to  indicate  which  one 
it  was.  Usually  termination  means  a  had  data  card  or  that  the  user  was 
trying  to  generate  points  beyond  the  capability  of  the  program. 

J.  COORDINATES 

The  cartesian  system  of  coordinates  x  and  y  are  u3ed  for  problems  in¬ 
volving  plane  geometry.  When  axisymmetric  problems  are  considered  the 
x,©-,z  system  of  polar  coordinates  is  used.  The  analysis  of  axisymmetric 
problems  is  conducted  in  the  x,z  plane. 

K.  TEMPERATURE  OUTPUT 

The  temperature  output  can  be  modified  in  two  ways.  The  number  of 
time  increments  printed  out  can  be  limited  by  using  the  control  parameter 
INTER.  If  the  program  is  directed  to  solve  for  the  temperature  every 
minute  for  IOC  minutes  and  INTER  is  equal  Co  25,  temperature  information 
will  only  be  printed  for  25,  50,  75,  and  100  minutes.  The  second  manner 
in  which  the  output  may  be  modified  is  at  the  end  of  each  time  sequence 
the  x,y  coordinates  of  each  node)  point  and  the  temperature  are  printed 
if  the  control  parameter  KPUNCH  equals  one. 

L.  INPUT  DATA  PREPARATION 


The  following  sections  describe  the  manner  in  which  the  input  data 
is  to  be  placed  on  data  cards.  The  data  is  prepared  in  groups  and  the 
final  data  deck  is  assembled  by  putting  the  data  groups  together  m  order 
of  their  listing  below. 
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1.  Control  Information 


(a)  Title  Card  (18A4):  72  arbitrary  spaces  tor  an  alphaaeric 

problem  label  punched  in  columns  1-72 

(b)  Units  Card  (12M):  Alphameric  unit  labels 

Columns  1-12  length 

13-24  m ass 

25- 36  time 

37-48  temperature 

Labels  should  be  left  justified. 

(c)  Control  Information  Card  (1115 ,  F10.0) 

Columns  1-5  NIHNP  number  of  nodal  points 

6-10  NOMEL  number  of  elements 

11-15  NIMCBC  number  of  convection 

boundary  conditions 
16-20  NUMMAT  number  of  materials 

21-25  NSEQ  number  of  time  sequences 

26- 30  KAT  type  of  geometry 

31-35  INTER  print  interval 

36-40  UNIT  units  system  (integer) 

41-45  KPUNCH  special  distribution 

46-50  temperature  conditions 

51-55  KEY  convection  conditions 

56-65  TO  initial  temperature 

"I”  formats  must  be  right  justified. 

2.  Material  Properties 

(a)  Material  Properties  Card  (5D10.3) 

Columns  1-10  XC0ND 
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11-20  YCOND  kyy 

21-30  XYCOND  kxy 

31-40  SPHT  fc 

41-50  HX  heat  generation 

One  card  for  each  material. 

3.  Time  Sequencing  Information 

(a)  Time  Sequencing  Card  (2L5,H0,D10.3): 

Columns  1-5  LT A3  temperature  dependent 

properties  (true  or  false) 
6-10  MTA3  time  dependent  boundary 

conditions  (true  or  false) 
11-2C  ITAG  number  of  time  steps 

21-21  TA3  time  increment 

One  card  for  each  time  sequence. 

4.  Nodal  Point  Information 

(a)  Nodal  Point  Card  (215, 3F 10.0) : 

Columns  1-5  N  nodal  point  number 

6-10  K0DE  temperature  boundary 

condition 

11-20  X  x  coordinate 

21-30  Y  y  coordinate 

31-40  T  initial  temperature 

One  card  for  each  nodal  point  not  generated  by  the  program. 

5.  Element  Information 

(a)  Element  Card  (615, D10. 3): 

Columns  1-5  N  element  number 

6-10  IX(N,1)  nodal  point  number 
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11-15 

IX (M, 2) 

nodal  point  number 

16-20 

IX  (N,  3) 

nodal  point  number 

21-25 

IX  (N,  4) 

nodal  point  number 

26-30 

IX(N,5) 

material  number 

31-40 

HTGEN 

element  heat  generation 

One  card  for  each  element  not  generated.  Nodal  points  are  ntaa- 

bered  clockwise  starting  in 

the  upper 

left  hand  corner. 

6. 

Convection  Boundarv  Conditions 

(a)  Convection  Boundary  Card  (215,010.3, F10.0) 

Columns  1-5 

IB 

nodal  point  number 

6-10 

JB 

nodal  point  number 

11-20 

HB 

convective  heat  transfer 

coefficient 

21-30 

TEMPB 

ambient  fluid  temperature 

One  card  for  each  convection  boundary  condition. 

7. 

(a)  Temperature  or  Flux  Boundary  Card  (I10,D10.3) 

Columns  1-10 

NN 

nodal  point  number 

11-20 

TQ 

temperature  or  flux 

One  card  for  each  tempe 

rature  or 

flux  boundary  point.  A  tempera 

ture  or  flux  boundary  card  for  the  last  nodal  point  oust  be  read  even 
if  it  isn’t  a  boundary  point. 

8.  Multiple  Problems 

The  above  cards  comprise  a  data  deck  for  one  problem.  Multiple 
problems  may  be  solved  by  putting  their  data  decks  together.  The  only 
limitation  on  this  procedure  is  computer  time.  It  is  important  that 
there  be  three  blank  cards  after  the  last  problem  data  set  in  the  date 
deck.  This  is  to  terminate  operation  of  the  program. 
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RECOMMENDATIONS 


Since  the  number  of  problem  types  that  can  be  solved  by  this  program 
is  limited  mostly  by  the  imagination  of  the  user  it  is  highly  recommended 
that  this  engineering  tool  should  be  brought  to  the  attention  of  those 
students  studying  in  the  area  of  unsteady  heat  conduction.  Only  through 
complete  familiarization  with  the  program  can  the  user  make  use  of  its 
great  potential  to  solve  problems  involving  complex  configurations  and 
boundary  conditions. 

Recommendations  for  further  study  in  this  area  are: 

1.  Expansion  of  the  finite  element  method  to  three  dimensional  con¬ 
figurations. 

2.  Revision  of  the  boundary  c;>  edition  segments  of  the  program  to 
make  use  of  random  access  disks  for  temporary  storage. 

3.  Construction  of  a  supplementary  program  that  will  punch  input 
data  cards  for  complex  time  dependent  boundary  conditions. 

4.  An  isothermal  contour  plotting  subroutine  capable  of  being  used 
in  conjunction  with  the  present  NFS  Computer  Center  subroutine  DRAW. 

5.  Use  of  the  program  in  conjunction  with  the  NPS  Computer  Center 
IBM  Optical  Display  Unit,  in  order  to  observe  transient  heat  conduction 
differences  after  minor  alterations  in  mesh  configuration. 

6.  Use  of  the  program  in  conjunction  with  existing  finite  element 
programs  in  the  area  of  stress  analysis  to  obtain  thermal  stresses. 
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APPENDIX  A 


In  this  section  some  solutions  to  different  configurations  will  be 
presented  as  examples  of  possible  problem  types  that  may  be  solved  with 
this  program.  Also  the  results  of  these  examples  are  compared  with 
analytic  solutions  of  the  same  configurations.  In.  order  to  avoid  con¬ 
fusion  the  same  material  nickel  (99.9%)  and  the  same  initial  temperature, 
100°F,  are  used  in  each  problem. 

The  material  properties  of  nickel  are: 

p  *»  556  lbm/ft^ 
c  -  .1065  BTU/lbm-°F 
«(T  «  .882  ft* /hr 

Example  1. 

This  is  a  problem  in  one  dimensional  heat  conduction.  It  could  also 
be  considered  a  problem  of  heat  conduction  in  a  send -infinite  plate  by 
assuming  that  lines  parallel  to  the  direction  of  heat  flux  are  adiabatic 
boundaries.  Because  of  this  assumption  we  only  need  one  row  of  elements 
across  the  plate  which  is  one  foot  wide  (see  figure  A-l),  The  top  and 
bottom  sides  of  the  elements. are  adiabatic  boundaries.  Also  the  end 
boundary  conditions  are  adiabatic  on  the  left  and  are  constant  0°F  on  the 
right.  The  nodal  points  are  numbered  from  left  to  right,  starting  with 

the  top  row.  The  elements  are  also  numbered  from  left  to  right.  The 
computer  solution  provides  answers  at  distinct  time  intervals  while  most 


Figure  A-j.  One  Dimensional  Heat  Conduction 
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analytic  solutions  are  given  in  terms  of  the  dimensionless  time  parameter 
Fourier  number,  Nj-q.  The  Fourier  number,  Npo,  equals  t/L  ,  where  **  t 
is  the  thermal  diffusivity  of  the  material.  The  comparisons  will  be  made 


at  various 

Fourier  numbers. 

The  analytic  solution  was 

obtained  from 

page  98  of 

reference  Q^J. 
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.4 

15 
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Example  2. 

This  example  has  the  same  mesh  configuration  and  boundary  conditions 


as  the  first  example.  The  initial  temperature  distribution  is  changed 
to  the  form  T  «  100  -  50x°F.  The  analytic  solution  was  obtained  from 
page  98  of  reference  DO- 
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Example  3. 

This  example  Is  a  solution  for  two-dimensional  heat  conduction.  The 
configuration  is  a  square  plate  vith  an  initial  temperature  of  100°F  and 
boundary  conditions  of  0°F  on  all  sides.  Using  symmetry  the  nser  need 
only  examine  one  fourth  of  the  plate,  i.e.  one  corner  (see  figure  A-2), 


Figure  A-2  Two  Dimensional  Heat  Conduction 


The  nodal  points  are  numbered  from  left  to  right  and  from  the  bottom  row 
to  the  top.  The  elements  are  numbered  in  the  same  manner.  In  this  case 
the  problem  consists  of  a  plate  with  two  sides  at  0°F  and  the  other  two 
insulated  for  an  adiabatic  boundary  condition.  The  analytic  solution  is 
obtained  by  using  the  solution  for  one  dimensional  neat  conduction  in  a 
semi -inf inite  plate  and  the  method  of  multiplicative  superposition. 
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The  analytic  aolution  war  obtained  from  page  118  of  reference  (jS}. 
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Example  4. 


In  this  problem  the  same  mesh  configuration  aa  the  preview  example 
is  used.  The  boundary  conditions  are  changed  so  that  the  sides  are  adia¬ 
batic,  the  top  is  at  a  constant  100°F.  The  Initial  temperature  is  100°F, 
By  taking  an  extremely  large  time  increment,  that  la  of  a  magnitude 
greater  than  the  maxitaum  number  of  significant  figures  available  in  the 
computer,  say  10^,  the  first  time  step  will  yield  the  steady  state 
aolution  to  the  problem.  This  is  made  possible  because  the  lime  depend¬ 
ent  contributions  to  the  coefficient  matrix  and  the  forcing  vector  are 
multiplied  by  the  Inverse  of  the  time  Increment  which  in  this  case  is 
essentially  saro  (10  ^ ).  Ths  solution  to  this  problem  is  a  linear 
tamparature  distribution  and  the  results  below  were  obtained  by  using  a 
time  increment  of  10^  minutes. 

Nodal  Finite 

Point  Element. 
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APPENDIX  B 


A 

AGAi/TR 

B 

CONDI 

CONDJ 

CONDJ 

K 

m 

?NTER 

hag 

IX 

JIMP 

KAT 

KEY 

KODK 

KPIJNCU 

LPRINT 

LTAG 

MB  AND 

MTA3 

MTYPE 

NSEQ 

NUMCBC 

NIMEL 

NUMMAT 


coefficient  matrix 
dummy  logicel  variable 
effective  forcing  vector 
k*  in  an  element 
ky  in  an  element 
kxy  in  an  element 

convection  boundary  layer  heat  transfer  coefficient 

heat  generation  in  a  material 

output-.  print  interval 

number  of  time  steps 

k  »  1  -  4  nodal  point  numbers 

k  «*  5  material  number 

temperature  and  flux  boundary  generation 
geometry  type 

convection  boundary  generation 

temperature  or  flux  boundary  parameter 

special  temperature  distribution  parameter 

print  parameter 

logical  variable 

half  band  width 

logical  variable 

material  number 

number  of  time  sequences 

number  of  convection  boundary  conditions 

number  of  elements 

number  of  materials 
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NUMNP 

SPHX 

T 

TJG 

TEMP 

TIME 

TITLE 

■MEAN 

TO 

UN 

UNIT 

X 

XCOND 

XMEAN 

XYCOND 

Y 

YCOND 

YMEAN 


number  of  nodal  points 

the  prodtict  fc 

nodal  temperatura 

time  increment 

ambient  fluid  temperature 

tin*  variable 

identification  title 

mean  temperature  in  an  element 

initial  temperature 

unit  labels 

unit  system  parameter 

nodal  point  x  coordinate 

k*  in  a  material 

mean  x  coordinate  of  an  element 

kxy  in  a  material 

nodal  point  y  coordinate 

ky  in  a  material 

mean  y  coordinate  of  an  element 
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APPENDIX  C 


* 

Title  and  Labels 

THESIS  CHECK  PROBLEM  HEAT  CONDUCTION  IN  A  RECTANGLE 
EEET  POUND  MASS  MINUTES  FAHRENHEIT 


*Control  Information 

36  25  0  1  1  C  I  6  1  0  0100,0 


*Material  Properties 

♦8* 66 00-01 ♦8,6600-0 1+0, 00 OD +00+5, 9300+01+0,0000+00 

*Time  Sequencing  Information 
FALSEFALSE  250+1.0000-01 


A 

Nodal  Point  Information 


l 

1C.0 

C.O 

100.0 

2 

10.1 

0.0 

100.0 

•a 

10.2 

0.0 

100.0 

4 

10.3 

0.0 

100.0 

5 

10.4 

0.0 

100.0 

6 

10.5 

0.0 

100.0 

7 

10.0 

0.1 

100.0 

12 

00.5 

0.1 

100.0 

13 

IC.O 

0.2 

100. 0 

18 

00.5 

0.2 

100,0 

19 

10.0 

0.3 

100. 0 

24 

00.5 

0.3 

100.0 

25 

10.0 

0.4 

100.0 

30 

00.5 

0.4 

100.0 

31 

10.0 

C.  5 

100.0 

36 

00.5 

0.5 

100.0 

Element  Information 


1 

7 

8 

2 

1 

1+0.0000+00 

5 

11 

12 

6 

5 

1+0. 0000+00 

6 

13 

1 A 

8 

7 

1+0.0000+00 

10 

17 

18 

12 

11 

1+0.0000+00 

11 

19 

20 

14 

13 

1+0. 0000+00 

15 

23 

24 

18 

17 

1+0.0000+00 

16 

25 

26 

20 

19 

1+0.0000+00 

20 

29 

30 

24 

23 

1+0.0000+00 

21 

31 

32 

26 

25 

1+0.0000+00 

25 

35 

36 

30 

29 

1+0.0000+00 

*Boundary  Conditions 

1+0, 0000+00 
2+0, OCCD+QO 
3+C. 0000+00 
4+0,0000+00 
5+O.OOCr+GO 
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